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It is shown both analytically and numerically that the suggested lattice Boltzmann model on 
rectangular grids leads to anisotropic dissipation of fluid momentum and thus it does not recover 
Navier-Stokes equations. Hence, it cannot be used for the simulation of hydrodynamics. 
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In a recent paper pj, a lattice Bhatnagar-Gross-Krook (LBGK) model is suggested for the simulation of the Navier- 
Stokes equation for incompressible flow. Setting itself apart from the standard LBGK equation 0] on square lattices, 
the LBGK model of Ref. |l| is formulated on generic rectangular lattices. This is surprising since earlier studies clearly 
indicated that it is impossible to formulate a single relaxation time LBGK model for the Navier-Stokes equation on 
simple rectangular lattices Q . 

However, a more careful study of the LBGK model in [l| detailed below reveals that the resulting hydrodynamic 
equations, in fact, lead to anisotropic viscous pressure tensor, thus disproving its main result. Indeed, following 
let us consider a rectangular two-dimensional lattice with the velocities e a = (e ax ,e ay ), a = 0, ...,8, where 
Sax 6 {—e x ,0,e x } and e ay £ {—e y ,0,e y } (the rectangular D2Q9 lattice, see Fig. 1 in [l|), and revisit the derivation 
of the Navier-Stokes equation from the LBGK model presented in Appendix B of [l|. Specifically, the computation 
of the non-equilibrium pressure tensor II^ 1 ^ (B23), (B24) amounts to evaluation of a function which appears in the 
left hand side of equation (B21), 

s a ie a je a kf a ^- (1) 
Let us denote the third-order moment at equilibrium in the above expression, 

Qijk = ^ ] e aiS-ajGakf a \ (2) 
a 

where is the equilibrium population / cq given by equation (11) of For the two-dimensional model there are 
four independent components of the tensor Q xxx , Qyyy, Q x ° yy and Q y ° xx . The two diagonal components, Q x a xx 
and Q y °J y , are independent of any particular form of the equilibrium because of the lattice constraint, e 3 ai = efe a i. 
Therefore, for the diagonal components we have Q xxx = pe%u x and Q yyy = pe y u y , irrespectively of the particular 
form of / oq , or, introducing the aspect ratio a = e x /e y , 

Qxxx = pe x e v au x , Q y ° y \ = pe x e y aT x u y (3) 
On the other hand, evaluating the remaining two off-diagonal components on the equilibrium (11), we find 

Qxyy = ^a-\ x , Q% = ^au y . (4) 

Note that equations ^ and (@| are exact. Thus, unless a = 1, tensor © is not isotropic, and hence the right 
hand side of (B21) as given in [l| is not correct. When the correct expressions ^ and ^ are used in the further steps 
of the derivation, we finally obtain, at variance with the isotropic Newtonian viscous stresses (B24), an anisotropic 
expression depending on the aspect ratio a, 
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where v is given by equation (B25) in Ref. If a = 1 (standard square grid), the viscous stresses ([5]) become 
isotropic, and v is interpreted as the kinematic viscosity. However, there is no such interpretation at a ^ 1, and ([5]) 
is not an isotropic viscous stress of the Navier-Stokes equation. Thus, the standard analysis of hydrodynamic limit 
reveals that the model [l[ does not recover Navier-Stokes equations at a / 1. 

Same consideration applies to the three-dimensional case of It is easy to see that a particular choice of the 
equilibrium is actually irrelevant because the anisotropy is already present in the diagonal components of the third- 
order moment which are independent of a particular choice of the equilibrium. 

To numerically probe the anisotropy in the model we study below simulations of the standard Taylor-Green 
(TG) vortex flow in two dimensions at Reynolds number Re = 10. The simulation with a standard square lattice is 
fully resolved and agrees well with the analytical solution on the grid with 32 x 32 nodes. However, when we refine 
the grid and increase the x resolution to 64 points while retaining the y resolution at 32 points (that is, e y = 2e x ), 
we clearly see the anisotropic behavior of viscosity. The velocity in the y direction decays much faster resulting 
in a deformed vortex which eventually corrupts the entire simulation. Fig. [T] also shows that the anisotropy may go 
unnoticed in the early stages of the simulation (see for example the snapshot at T = 500 in Fig. [T]) . In the simulation, 
Reynolds number was defined as Re = UqN x /v 1 where Uq is the amplitude of the velocity of TG vortex, N x is the 
number of nodes in the x direction and v is the parameter in ([5]) related to the LBGK relaxation time r according to 
equation (B25) of [l|: v = e x e y (2r — l)/6. 

We note in passing that at a = 1 the equilibrium (11) in [l| is different from the standard equilibrium used in 
LBGK models, in particular, equilibrium populations of the diagonal velocities become equal to zero at vanishing flow 
velocity. This probably explains why model [l[ is numerically less stable than the standard LBGK model on a square 
lattice. The simulation of Taylor-Green vortex flow on a grid of 32 x 32 remains stable for up to Re = 10000 for the 
standard lattice Boltzmann scheme while it fails below Re = 150 with the scheme pj. 

Finally, a comment on the simulations presented in [l[ is in order. For some flow configurations (unidirectional 
flows) the above anisotropy may be not very visible in the velocity plots. All simulations with rectangular lattices in 
[l[ concern predominantly unidirectional flows without vortex structures. Such flows do not provide a benchmark for 
an anisotropic model of |l| as they exclude rotation of the flow. Even the simplest vortex simulation reveals unphysical 
anisotropy of the model, as was shown above. 

To conclude, the model suggested in [l[ does not recover the Navier-Stokes equations for non-square rectangular 
lattices and cannot be used for simulation of hydrodynamics. 
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(b) T = 500 
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(c) T = 1100 (d) T = 1500 

FIG. 1: Taylor-Green vortex simulation with model [l[ at Reynolds number Re = 10, initial velocity Uq = 0.05 and a 
grid of N x = 64, N y = 32 (e y = 2e x ), v — Uq x N x /Re = 0.32. Velocity vector plots (with uniform vector length) at 
various times show the accumulation of anisotropy: The vortex is visibly corrupted at T = 1100 and is fully 

destroyed at T = 1500. 



